###############
### Setup R ###
###############

sink("results.txt")

### Clear space
rm(list = ls())

### Load packages
library(rio)
library(data.table)
library(scales)

### Print system info
sessionInfo()

##################
### Setup Data ###
##################

### Load data
dat <- rio::import("copecrabtreefariss.csv")
dim(dat)

#################################################
### Create Main Manuscript Tables and Figures ###
#################################################

print("Table 3.2 results")
### Create Table 3.2
#### First row
tort <- na.omit(data.frame(dat$year, dat$ccode, dat$CTRY, dat$TORT, dat$v2cltort))
dim(tort)[1] ##### N
na.tort <- dim(dat)[1] - dim(tort)[1]
na.tort ##### Dropped 
cor.tort <- cor(tort$dat.TORT, tort$dat.v2cltort, method = c("spearman"))
round(cor.tort, 3) ##### Rho
#### Second row
kill <- na.omit(data.frame(dat$year, dat$ccode, dat$CTRY, dat$KILL, dat$v2clkill))
dim(kill)[1] ##### N
na.kill <- dim(dat)[1] - dim(kill)[1]
na.kill ##### Dropped
cor.kill <- cor(kill$dat.KILL, kill$dat.v2clkill, method = c("spearman"))
round(cor.kill, 3) ##### Rho
#### Third row
assn <- na.omit(data.frame(dat$year, dat$ccode, dat$CTRY, dat$ASSN, dat$v2x_frassoc_thick))
dim(assn)[1] ##### N
na.assn <- dim(dat)[1] - dim(assn)[1]
na.assn ##### Dropped
cor.assn <- cor(assn$dat.ASSN, assn$dat.v2x_frassoc_thick, method = c("spearman"))
round(cor.assn, 3)
#### Fourth row
dmov <- na.omit(data.frame(dat$year, dat$ccode, dat$CTRY, dat$DOMMOV, dat$v2xcl_dmove))
dim(dmov)[1] ##### N
na.dmov <- dim(dat)[1] - dim(dmov)[1]
na.dmov ##### Dropped
cor.dm <- cor(dmov$dat.DOMMOV, dmov$dat.v2xcl_dmove, method = c("spearman"))
round(cor.dm, 3)
#### Fifth row
fmov <- na.omit(data.frame(dat$year, dat$ccode, dat$CTRY, dat$FORMOV, dat$v2clfmove))
dim(fmov)[1] ##### N
na.fmov <- dim(dat)[1] - dim(fmov)[1]
na.fmov ##### Dropped
cor.fm <- cor(fmov$dat.FORMOV, fmov$dat.v2clfmove, method = c("spearman"))
round(cor.fm, 3)
#### Sixth row
rel <- na.omit(data.frame(dat$year, dat$ccode, dat$CTRY, dat$NEW_RELFRE, dat$v2clrelig))
dim(rel)[1] ##### N
na.rel <- dim(dat)[1] - dim(rel)[1]
na.rel ##### Dropped
cor.rel <- cor(rel$dat.NEW_RELFRE, rel$dat.v2clrelig, method = c("spearman"))
round(cor.rel, 3)
##### Average reported in text 
round((0.471 + 0.526 + 0.741 + 0.476 + 0.585 + 0.575)/6, 3)

### Create Figure 1
#### Setup data
kill <- data.table(kill)
dat.kill.year.cor <- kill[, .(mCor = cor(dat.KILL, dat.v2clkill, method = c("spearman"))), by = dat.year]
dat.kill.year.cor
tort <- data.table(tort)
dat.tort.year.cor <- tort[, .(mCor = cor(dat.TORT, dat.v2cltort, method = c("spearman"))), by = dat.year]
dat.tort.year.cor
assn <- data.table(assn)
dat.assn.year.cor <- assn[, .(mCor = cor(dat.ASSN, dat.v2x_frassoc_thick, method = c("spearman"))), by = dat.year]
dat.assn.year.cor
rel <- data.table(rel)
dat.rel.year.cor <- rel[, .(mCor = cor(dat.NEW_RELFRE, dat.v2clrelig, method = c("spearman"))), by = dat.year]
dat.rel.year.cor
dmov <- data.table(dmov)
dat.dmov.year.cor <- dmov[, .(mCor = cor(dat.DOMMOV, dat.v2xcl_dmove, method = c("spearman"))), by = dat.year]
dat.dmov.year.cor
fmov <- data.table(fmov)
dat.fmov.year.cor <- fmov[, .(mCor = cor(dat.FORMOV, dat.v2clfmove, method = c("spearman"))), by = dat.year]
dat.fmov.year.cor

#### Create panel of plots
png("yearcor.png", width = 1000, height = 1100)
par(mfrow = c(3, 2), oma = c(2, 2, 2, 2), mar = c(5.5, 4.5, 2.5, 1))
plot(dat.kill.year.cor$dat.year, dat.kill.year.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(0.4, 1),
     ylab = "Correlation Coefficient", xlab = "Year", 
     main = "Extrajudicial Killing", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, 
     cex.main = 2, cex.sub = 1.25)
abline(lm(dat.kill.year.cor$mCor ~ dat.kill.year.cor$dat.year), col = "#666666") 
lines(lowess(dat.kill.year.cor$mCor ~ dat.kill.year.cor$dat.year), col = "#333333", lwd = 2)
plot(dat.tort.year.cor$dat.year, dat.tort.year.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(0.4, 1),
     ylab = "Correlation Coefficient", xlab = "Year", 
     main = "Torture", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, 
     cex.main = 2, cex.sub = 1.25)
abline(lm(dat.tort.year.cor$mCor ~ dat.tort.year.cor$dat.year), col = "#666666") 
lines(lowess(dat.tort.year.cor$mCor ~ dat.tort.year.cor$dat.year), col = "#333333", lwd = 2)
plot(dat.assn.year.cor$dat.year, dat.assn.year.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(0.4, 1),
     ylab = "Correlation Coefficient", xlab = "Year", 
     main = "Association", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, 
     cex.main = 2, cex.sub = 1.25)
abline(lm(dat.assn.year.cor$mCor ~ dat.assn.year.cor$dat.year), col = "#666666") 
lines(lowess(dat.assn.year.cor$mCor ~ dat.assn.year.cor$dat.year), col = "#333333", lwd = 2)
plot(dat.rel.year.cor$dat.year, dat.rel.year.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(0.4, 1),
     ylab = "Correlation Coefficient", xlab = "Year", 
     main = "Religious Freedom", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, 
     cex.main = 2, cex.sub = 1.25)
abline(lm(dat.rel.year.cor$mCor ~ dat.rel.year.cor$dat.year), col = "#666666") 
lines(lowess(dat.rel.year.cor$mCor ~ dat.rel.year.cor$dat.year), col = "#333333", lwd = 2)
plot(dat.dmov.year.cor$dat.year, dat.dmov.year.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(0.4, 1),
     ylab = "Correlation Coefficient", xlab = "Year", 
     main = "Domestic Movement", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, 
     cex.main = 2, cex.sub = 1.25)
abline(lm(dat.dmov.year.cor$mCor ~ dat.dmov.year.cor$dat.year), col = "#666666") 
lines(lowess(dat.dmov.year.cor$mCor ~ dat.dmov.year.cor$dat.year), col = "#333333", lwd = 2)
plot(dat.fmov.year.cor$dat.year, dat.fmov.year.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(0.4, 1),
     ylab = "Correlation Coefficient", xlab = "Year", 
     main = "Foreign Movement", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, 
     cex.main = 2, cex.sub = 1.25)
abline(lm(dat.fmov.year.cor$mCor ~ dat.fmov.year.cor$dat.year), col = "#666666") 
lines(lowess(dat.fmov.year.cor$mCor ~ dat.fmov.year.cor$dat.year), col = "#333333", lwd = 2)
dev.off()

### Create Table 3.3
#### Setup data, building on data used to create Figure 1
dat.kill.country.cor <- kill[, .(mCor = cor(dat.KILL, dat.v2clkill, method = c("spearman"))), by = dat.CTRY]
dat.kill.country.cor
dat.tort.country.cor <- tort[, .(mCor = cor(dat.TORT, dat.v2cltort, method = c("spearman"))), by = dat.CTRY]
dat.tort.country.cor
dat.assn.country.cor <- assn[, .(mCor = cor(dat.ASSN, dat.v2x_frassoc_thick, method = c("spearman"))), by = dat.CTRY]
dat.assn.country.cor
dat.dmov.country.cor <- dmov[, .(mCor = cor(dat.DOMMOV, dat.v2xcl_dmove, method = c("spearman"))), by = dat.CTRY]
dat.dmov.country.cor
dat.fmov.country.cor <- fmov[, .(mCor = cor(dat.FORMOV, dat.v2clfmove, method = c("spearman"))), by = dat.CTRY]
dat.fmov.country.cor
dat.rel.country.cor <- rel[, .(mCor = cor(dat.NEW_RELFRE, dat.v2clrelig, method = c("spearman"))), by = dat.CTRY]
dat.rel.country.cor

print("Table 3.3 results")
#### First row
dim(dat.kill.country.cor[dat.kill.country.cor$mCor < 0, ])[1] ##### Neg. N
dim(dat.kill.country.cor)[1] ##### total N
round(dim(dat.kill.country.cor[dat.kill.country.cor$mCor < 0, ])[1] / dim(dat.kill.country.cor)[1], 2)*100 ##### % Neg.
#### Second row
dim(dat.tort.country.cor[dat.tort.country.cor$mCor < 0, ])[1] ##### Neg. N
dim(dat.tort.country.cor)[1] ##### total N
round(dim(dat.tort.country.cor[dat.tort.country.cor$mCor < 0, ])[1] / dim(dat.tort.country.cor)[1], 2)*100 ##### % Neg.
#### Third row
dim(dat.assn.country.cor[dat.assn.country.cor$mCor < 0, ])[1] ##### Neg. N
dim(dat.assn.country.cor)[1] ##### total N
round(dim(dat.assn.country.cor[dat.assn.country.cor$mCor < 0, ])[1] / dim(dat.assn.country.cor)[1], 2)*100 ##### % Neg.
#### Fourth row
dim(dat.dmov.country.cor[dat.dmov.country.cor$mCor < 0, ])[1] ##### Neg. N
dim(dat.dmov.country.cor)[1] ##### total N
round(dim(dat.dmov.country.cor[dat.dmov.country.cor$mCor < 0, ])[1] / dim(dat.dmov.country.cor)[1], 2)*100 ##### % Neg.
#### Fifth row
dim(dat.fmov.country.cor[dat.fmov.country.cor$mCor < 0, ])[1] ##### Neg. N
dim(dat.fmov.country.cor)[1] ##### total N
round(dim(dat.fmov.country.cor[dat.fmov.country.cor$mCor < 0, ])[1] / dim(dat.fmov.country.cor)[1], 2)*100 ##### % Neg.
#### Sixth row
dim(dat.rel.country.cor[dat.rel.country.cor$mCor < 0, ])[1] ##### Neg. N
dim(dat.rel.country.cor)[1] ##### total N
round(dim(dat.rel.country.cor[dat.rel.country.cor$mCor < 0, ])[1] / dim(dat.rel.country.cor)[1], 2)*100 ##### % Neg.
##### Empowerment rights average reported in text 
round((20 + 29 + 23 + 29)/4, 0)

### Create Figure 2
png("countrycor.png", width = 1000, height = 1100)
par(mfrow = c(3, 2), oma = c(2, 2, 2, 2), mar = c(5.5, 4.5, 2.5, 1))
plot(as.numeric(dat.kill.country.cor$dat.CTRY), dat.kill.country.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(-1, 1),
     ylab = "Correlation Coefficient", xlab = "Country Name", 
     main = "Extrajudicial Killing", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, xaxt = 'n',
     cex.main = 2, cex.sub = 1.25, frame.plot = F)
axis(1, at = as.numeric(dat.kill.country.cor$dat.CTRY), labels = dat.kill.country.cor$dat.CTRY)
abline(h = 0, col = "#999999", lty = 3, lwd = 3)
plot(as.numeric(dat.tort.country.cor$dat.CTRY), dat.tort.country.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(-1, 1),
     ylab = "Correlation Coefficient", xlab = "Country Name", 
     main = "Torture", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, xaxt = 'n', 
     cex.main = 2, cex.sub = 1.25, frame.plot = F)
axis(1, at = as.numeric(dat.tort.country.cor$dat.CTRY), labels = dat.tort.country.cor$dat.CTRY)
abline(h = 0, col = "#999999", lty = 3, lwd = 3)
plot(as.numeric(dat.assn.country.cor$dat.CTRY), dat.assn.country.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(-1, 1),
     ylab = "Correlation Coefficient", xlab = "Country Name", 
     main = "Association", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, xaxt = 'n',
     cex.main = 2, cex.sub = 1.25, frame.plot = F)
axis(1, at = as.numeric(dat.assn.country.cor$dat.CTRY), labels = dat.assn.country.cor$dat.CTRY)
abline(h = 0, col = "#999999", lty = 3, lwd = 3)
plot(as.numeric(dat.rel.country.cor$dat.CTRY), dat.rel.country.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(-1, 1),
     ylab = "Correlation Coefficient", xlab = "Country Name", 
     main = "Religious Freedom", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, xaxt = 'n',
     cex.main = 2, cex.sub = 1.25, frame.plot = F)
axis(1, at = as.numeric(dat.rel.country.cor$dat.CTRY), labels = dat.rel.country.cor$dat.CTRY)
abline(h = 0, col = "#999999", lty = 3, lwd = 3)
plot(as.numeric(dat.dmov.country.cor$dat.CTRY), dat.dmov.country.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(-1, 1),
     ylab = "Correlation Coefficient", xlab = "Country Name", 
     main = "Domestic Movement", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, xaxt = 'n',
     cex.main = 2, cex.sub = 1.25, frame.plot = F)
axis(1, at = as.numeric(dat.dmov.country.cor$dat.CTRY), labels = dat.dmov.country.cor$dat.CTRY)
abline(h = 0, col = "#999999", lty = 3, lwd = 3)
plot(as.numeric(dat.fmov.country.cor$dat.CTRY), dat.fmov.country.cor$mCor, 
     bty = 'n', col = "black", 
     pch = 19, ylim = c(-1, 1),
     ylab = "Correlation Coefficient", xlab = "Country Name", 
     main = "Foreign Movement", tck = -0.01,
     cex.lab = 1.75, cex.axis = 1.5, xaxt = 'n',
     cex.main = 2, cex.sub = 1.25, frame.plot = F)
axis(1, at = as.numeric(dat.fmov.country.cor$dat.CTRY), labels = dat.fmov.country.cor$dat.CTRY)
abline(h = 0, col = "#999999", lty = 3, lwd = 3)
dev.off()

### Create Table 3.4
#### Second column
dat.kill.country.cor[dat.kill.country.cor$dat.CTRY == "United States of America", ]
dat.tort.country.cor[dat.tort.country.cor$dat.CTRY == "United States of America", ]
dat.assn.country.cor[dat.assn.country.cor$dat.CTRY == "United States of America", ]
dat.dmov.country.cor[dat.dmov.country.cor$dat.CTRY == "United States of America", ]
dat.fmov.country.cor[dat.fmov.country.cor$dat.CTRY == "United States of America", ]
dat.rel.country.cor[dat.rel.country.cor$dat.CTRY == "United States of America", ]
#### Third column
dat.kill.country.cor[dat.kill.country.cor$dat.CTRY == "China", ]
dat.tort.country.cor[dat.tort.country.cor$dat.CTRY == "China", ]
dat.assn.country.cor[dat.assn.country.cor$dat.CTRY == "China", ]
dat.dmov.country.cor[dat.dmov.country.cor$dat.CTRY == "China", ]
dat.fmov.country.cor[dat.fmov.country.cor$dat.CTRY == "China", ]
dat.rel.country.cor[dat.rel.country.cor$dat.CTRY == "China", ]
#### Fourth column
dat.kill.country.cor[dat.kill.country.cor$dat.CTRY == "Germany", ]
dat.tort.country.cor[dat.tort.country.cor$dat.CTRY == "Germany", ]
dat.assn.country.cor[dat.assn.country.cor$dat.CTRY == "Germany", ]
dat.dmov.country.cor[dat.dmov.country.cor$dat.CTRY == "Germany", ]
dat.fmov.country.cor[dat.fmov.country.cor$dat.CTRY == "Germany", ]
dat.rel.country.cor[dat.rel.country.cor$dat.CTRY == "Germany", ]
#### Fifth column
dat.kill.country.cor[dat.kill.country.cor$dat.CTRY == "Russia", ]
dat.tort.country.cor[dat.tort.country.cor$dat.CTRY == "Russia", ]
dat.assn.country.cor[dat.assn.country.cor$dat.CTRY == "Russia", ]
dat.dmov.country.cor[dat.dmov.country.cor$dat.CTRY == "Russia", ]
dat.fmov.country.cor[dat.fmov.country.cor$dat.CTRY == "Russia", ]
dat.rel.country.cor[dat.rel.country.cor$dat.CTRY == "Russia", ]
#### Sixth column
dat.kill.country.cor[dat.kill.country.cor$dat.CTRY == "South Africa", ]
dat.tort.country.cor[dat.tort.country.cor$dat.CTRY == "South Africa", ]
dat.assn.country.cor[dat.assn.country.cor$dat.CTRY == "South Africa", ]
dat.dmov.country.cor[dat.dmov.country.cor$dat.CTRY == "South Africa", ]
dat.fmov.country.cor[dat.fmov.country.cor$dat.CTRY == "South Africa", ]
dat.rel.country.cor[dat.rel.country.cor$dat.CTRY == "South Africa", ]
#### Seventh column
dat.kill.country.cor[dat.kill.country.cor$dat.CTRY == "Iran", ]
dat.tort.country.cor[dat.tort.country.cor$dat.CTRY == "Iran", ]
dat.assn.country.cor[dat.assn.country.cor$dat.CTRY == "Iran", ]
dat.dmov.country.cor[dat.dmov.country.cor$dat.CTRY == "Iran", ]
dat.fmov.country.cor[dat.fmov.country.cor$dat.CTRY == "Iran", ]
dat.rel.country.cor[dat.rel.country.cor$dat.CTRY == "Iran", ]

#################################################
### Create Online Appendix Tables and Figures ###
#################################################

### Create Online Appendix A
dim(dat.kill.country.cor[!is.na(dat.kill.country.cor$mCor), ])[1]
dim(dat.tort.country.cor[!is.na(dat.tort.country.cor$mCor), ])[1]
dim(dat.assn.country.cor[!is.na(dat.assn.country.cor$mCor), ])[1]
dim(dat.dmov.country.cor[!is.na(dat.dmov.country.cor$mCor), ])[1]
dim(dat.fmov.country.cor[!is.na(dat.fmov.country.cor$mCor), ])[1]
dim(dat.rel.country.cor[!is.na(dat.rel.country.cor$mCor), ])[1]
#### Second column
dim(dat.kill.country.cor)[1] - dim(dat.kill.country.cor[!is.na(dat.kill.country.cor$mCor), ])[1]
dim(dat.tort.country.cor)[1] - dim(dat.tort.country.cor[!is.na(dat.tort.country.cor$mCor), ])[1]
dim(dat.assn.country.cor)[1] - dim(dat.assn.country.cor[!is.na(dat.assn.country.cor$mCor), ])[1]
dim(dat.dmov.country.cor)[1] - dim(dat.dmov.country.cor[!is.na(dat.dmov.country.cor$mCor), ])[1]
dim(dat.fmov.country.cor)[1] - dim(dat.fmov.country.cor[!is.na(dat.fmov.country.cor$mCor), ])[1]
dim(dat.assn.country.cor)[1] - dim(dat.rel.country.cor[!is.na(dat.rel.country.cor$mCor), ])[1]

### Create Online Appendix B
#### Prepare and check data
dat.kill.country.cor <- dat.kill.country.cor[order(dat.kill.country.cor$mCor), ]
dat.kill.country.cor <- dat.kill.country.cor[!is.na(dat.kill.country.cor$mCor), ]
length(unique(dat.kill.country.cor$dat.CTRY))
(length(dat.kill.country.cor$mCor) - length(unique(dat.kill.country.cor$dat.CTRY))) == 0
#### Make plot
png(height = 2000, width = 1300, pointsize = 20, file = "kill_ctry.png")
dotchart(dat.kill.country.cor$mCor, labels = dat.kill.country.cor$dat.CTRY, 
         pch = 19, lcolor = "white", main = "Extrajudicial Killing", 
         xlab = "Correlation Coefficient", cex = 0.7, cex.lab = 2, 
         xlim = c(-1, 1), sub = paste0("Number of Countries: ", length(unique(dat.kill.country.cor$dat.CTRY))))
abline(v = 0, col = alpha("#999999", .5), lty = 3, lwd = 3)
dev.off()

### Create Online Appendix C
#### Prepare and check data
dat.tort.country.cor <- dat.tort.country.cor[order(dat.tort.country.cor$mCor), ]
dat.tort.country.cor <- dat.tort.country.cor[!is.na(dat.tort.country.cor$mCor), ]
length(unique(dat.tort.country.cor$dat.CTRY))
(length(dat.tort.country.cor$mCor) - length(unique(dat.tort.country.cor$dat.CTRY))) == 0
#### Make plot
png(height = 2000, width = 1300, pointsize = 20, file = "torture_ctry.png")
dotchart(dat.tort.country.cor$mCor, labels = dat.tort.country.cor$dat.CTRY, pch = 19, 
         lcolor = "white", main = "Torture", xlab = "Correlation Coefficient", cex = 0.7, cex.lab = 2, xlim = c(-1, 1), 
         sub = paste0("Number of Countries: ", length(unique(dat.tort.country.cor$dat.CTRY))))
abline(v = 0, col = alpha("#999999", .5), lty = 3, lwd = 3)
dev.off()

### Create Online Appendix D
#### Prepare and check data
dat.assn.country.cor <- dat.assn.country.cor[order(dat.assn.country.cor$mCor), ]
dat.assn.country.cor <- dat.assn.country.cor[!is.na(dat.assn.country.cor$mCor), ]
length(unique(dat.assn.country.cor$dat.CTRY))
(length(dat.assn.country.cor$mCor) - length(unique(dat.assn.country.cor$dat.CTRY))) == 0
#### Make plot
png(height = 2000, width = 1300, pointsize = 20, file = "assn_ctry.png")
dotchart(dat.assn.country.cor$mCor, labels = dat.assn.country.cor$dat.CTRY, 
         pch = 19, lcolor = "white", main = "Association", 
         xlab = "Correlation Coefficient", cex = 0.7, cex.lab = 2, 
         xlim = c(-1, 1), sub = paste0("Number of Countries: ", length(unique(dat.assn.country.cor$dat.CTRY))))
abline(v = 0, col = alpha("#999999", .5), lty = 3, lwd = 3)
dev.off()

### Create Online Appendix E
#### Prepare and check data
dat.dmov.country.cor <- dat.dmov.country.cor[order(dat.dmov.country.cor$mCor), ]
dat.dmov.country.cor <- dat.dmov.country.cor[!is.na(dat.dmov.country.cor$mCor), ]
length(unique(dat.dmov.country.cor$dat.CTRY))
(length(dat.dmov.country.cor$mCor) - length(unique(dat.dmov.country.cor$dat.CTRY))) == 0
#### Make plot
png(height = 2000, width = 1300, pointsize = 20, file = "dm_ctry.png")
dotchart(dat.dmov.country.cor$mCor, labels = dat.dmov.country.cor$dat.CTRY, 
         pch = 19, lcolor = "white", main = "Domestic Movement", 
         xlab = "Correlation Coefficient", cex = 0.7, cex.lab = 2, 
         xlim = c(-1, 1), sub = paste0("Number of Countries: ", length(unique(dat.dmov.country.cor$dat.CTRY))))
abline(v = 0, col = alpha("#999999", .5), lty = 3, lwd = 3)
dev.off()

### Create Online Appendix F
#### Prepare and check data
dat.fmov.country.cor <- dat.fmov.country.cor[order(dat.fmov.country.cor$mCor), ]
dat.fmov.country.cor <- dat.fmov.country.cor[!is.na(dat.fmov.country.cor$mCor), ]
length(unique(dat.fmov.country.cor$dat.CTRY))
(length(dat.fmov.country.cor$mCor) - length(unique(dat.fmov.country.cor$dat.CTRY))) == 0
#### Make plot
png(height = 2000, width = 1300, pointsize = 20, file = "fm_ctry.png")
dotchart(dat.fmov.country.cor$mCor, labels = dat.fmov.country.cor$dat.CTRY, 
         pch = 19, lcolor = "white", main = "Foreign Movement", 
         xlab = "Correlation Coefficient", cex = 0.7, cex.lab = 2, 
         xlim = c(-1, 1), sub = paste0("Number of Countries: ", length(unique(dat.fmov.country.cor$dat.CTRY))))
abline(v = 0, col = alpha("#999999", .5), lty = 3, lwd = 3)
dev.off()

### Create Online Appendix G
#### Prepare and check data
dat.rel.country.cor <- dat.rel.country.cor[order(dat.rel.country.cor$mCor), ]
dat.rel.country.cor <- dat.rel.country.cor[!is.na(dat.rel.country.cor$mCor), ]
length(unique(dat.rel.country.cor$dat.CTRY))
(length(dat.rel.country.cor$mCor) - length(unique(dat.rel.country.cor$dat.CTRY))) == 0
#### Make plot
png(height = 2000, width = 1300, pointsize = 20, file = "religion_ctry.png")
dotchart(dat.rel.country.cor$mCor, labels = dat.rel.country.cor$dat.CTRY, 
         pch = 19, lcolor = "white", main = "Religious Freedom", 
         xlab = "Correlation Coefficient", cex = 0.7, cex.lab = 2, 
         xlim = c(-1, 1), sub = paste0("Number of Countries: ", length(unique(dat.rel.country.cor$dat.CTRY))))
abline(v = 0, col = alpha("#999999", .5), lty = 3, lwd = 3)
dev.off()

sink()